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ABSTRACT 

This paper presents an analytical formula that closely approximates the 
fully nonlinear power spectrum of matter fluctuations for redshift z ~ 5 to 
over a wide range of cosmologically interesting flat models with varying matter 
density Q m and neutrino fraction Q u . The functional form is motivated by 
analytical solutions in asymptotic regimes, but in order to obtain accurate 
approximations, the coefficients are calculated from fits to the nonlinear power 
spectrum computed from numerical simulations of four cosmological models. 
The transformation from the linear to nonlinear power spectrum depends on fl m , 
Q u , and time. A simple scaling rule is introduced, which greatly simplifies the 
construction of the functional form and allows the formula to depend directly 
on the rms linear mass fluctuation cr 8 instead of on an effective spectral index as 
in earlier work. 

Subject headings: cosmology : theory - dark matter - elementary particles - 
large-scale structure of universe - methods: analytical 
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1. Introduction 

The power spectrum of matter fluctuations P{k) provides the most fundamental 
statistical measure of gravitational clustering. When the amplitude of the density 
fluctuations is small, the power spectrum can be calculated easily from the linear 
perturbation theory of gravitational collapse. In the nonlinear regime, however, 
determination of the fully evolved power spectrum for a given cosmological model requires 
numerical simulations. Since nearly all observable astronomical systems have experienced 
some nonlinear collapse, it would provide much physical insight and practical convenience 
to devise a general analytical approximation (based on simulation results) for the nonlinear 
power spectrum for a wide range of cosmologically interesting models. 

By fitting to iV-body results, Hamilton et al. (1991) studied scale-free models with a 
power-law spectrum and presented a simple analytical formula that relates the spatially 
averaged nonlinear and linear two-point correlation function, £ n i(r) and £i(r ), where r 
is related to its pre-collapsed linear scale r by r = r (1 + £ n i) 1,/3 - This transformation 
then appeared to be magically independent of the spectral index n assumed in the model. 
Further tests against numerical simulations, however, found significant errors when the 
Hamilton et al. function was applied to models with n < — 1 (Jain, Mo, & White 1995; 
Padmanabhan 1996). Jain et al. (1995) instead proposed n-dependent formulas to relate £ 
and P(k) in the linear and nonlinear regimes in both scale-free models and the standard 
CDM model. For the more realistic CDM model, for which the spectral index changes 
from the primordial value n m 1 on large scales to nearly —3 on small scales, they used 
an effective index given by n e s = din P(k)/d\n k\k c , where k c is the scale at which the 
rms mass fluctuation o is unity. The index n e g therefore reflects the slope of the power 
spectrum at the length scale where nonlinearity becomes important. Peacock and Dodds 
(1996) extended this work to allow for a low f2 m and a non-zero cosmological constant. 
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No previous work has investigated in detail the subject of linear to nonlinear mapping 
in cold+hot dark matter (C+HDM) models that assume massive neutrinos are a component 
of the dark matter. This is perhaps because the physics in C+HDM models is generally 
more complicated than in CDM or LCDM (CDM with a cosmological constant) models 
due to the additional length scale associated with the free streaming of the neutrinos (Ma 
1996). Nevertheless, massive neutrinos remain a prime dark matter candidate, and the 
recent evidence for neutrino masses from the Super-Kamiokande experiment has made 
this possibility particularly intriguing (Fukuda et al. 1998). Although neither Jain et al. 
or Peacock & Dodds has tested these models, one may surmise that their formulas can 
be naturally extended to C+HDM models as long as the spectral index in the formula is 
calculated from the C+HDM power spectrum. This unfortunately does not work. Both 
fitting functions underestimate the nonlinear density variance^] A n i = Awk 3 P n \(k), at 
k^,2h Mpc -1 in C+HDM models, and the errors in Peacock & Dodds , for example, reach 
~ 50% at k ~ 10 h Mpc -1 . The linear to nonlinear transformation is therefore regulated 
by more than simply n e g. (Smith et al. (1998) recently reported agreement between 
the Peacock-Dodds formula and results from two C+HDM simulations. Their simulation 
resolution of ~ 0.3 h~ x Mpc, however, limited their test to only the mildly nonlinear regime, 
and could not probe the nonlinear regime where the large discrepancies reside.) 

This Letter differs from previous work in two ways. First, the simple analytical formula 
presented here closely approximates the fully nonlinear power spectrum of mass fluctuations 
at z & 5 in the previously unexplored C+HDM models as well as LCDM models with 
varying Q m . Numerical simulations of four COBE-normalized flat C+HDM and flat LCDM 
models are performed to calibrate the coefficients in the analytical formula. Second, the 

x The notations A and P(k) here are the same as in Jain et al., and are equivalent to A 2 
and P{k)/(2n) 3 in Peacock & Dodds . 
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formula introduced here depends directly on <j 8 (the rms linear mass fluctuation on 8 Yr 1 
Mpc scale) instead of a spectral index as in previous work. This is achieved by recognizing 
a scaling rule (see § 3), which also greatly simplifies the construction of the analytical 
formula. This work also extends into the nonlinear regime a previous investigation of the 
effects of neutrino free-streaming on the linear C+HDM power spectrum (Ma 1996). 



2. Input Linear Power Spectrum 

For a wide range of CDM and LCDM models that assume neutrinos are massless, a 
good approximation to the linear power spectrum is given by 

p(kn o _ m _ Ak-[D{a)/D,f [\n{l + a iq )/a iq f 

where k is the wavenumber in units of Mpc -1 , q = k/T h, T is a shape parameter, 
and a± = 2.34, a 2 = 3.89, 0:3 = 16.1,0:4 = 5.46, and 0:5 = 6.71 (Bardeen et 
al. 1986). The shape parameter T characterizes the dependence on cosmological 
parameters and is well approximated by T = Q m h/ exp[f2 b (l + l/f2 m )] (Efstathiou 
et al. 1992; Sugiyama 1995; see also Bunn & White 1997). The function D(a) is 
the linear growth factor, whose present value is D = D(a = 1), and it can be 
expressed as D(a) = ag, where the relative growth factor g is well approximated by 
g(n m (a), Q A (a)) = 2.5 n ni (a)[n m (a) 4 / 7 - Q A (a) + (1 + ft m (a)/2) (1 + fi A (a)/70)]- 1 (Lahav 
et al. 1991; Carroll et al. 1992). In LCDM models, g ~ 1 until the universe becomes 
A-dominated at 1 + z w ^m 1 ^ 3 ; the value of g then decreases with increasing a. The 
normalization factor A can be chosen by fixing the value of <t 8 ; if instead the COBE 
normalization is desired, it is A = 5 2 H (c/ 'H ) n+3 / '(47r) , where (for flat models) 5h = 
1.94 x 1(T 5 fi m - 785 - - 05 ln nm exp(-0.95n - 0.169n 2 ) with n = n - 1 (Bunn & White 1997). 

The linear power spectra for the C+HDM models require additional treatment since 
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the effect of massive neutrinos on the shape of the power spectrum is both time and scale 
dependent. It is found that by introducing a second shape parameter, T u = a}l 2 VL v h 2 , to 
characterize the neutrino free-streaming distance, one can obtain a good approximation 
to the linear power spectra (density averaged over the cold and hot components) in flat 
C+HDM models at z <, 5 when neutrinos are adequately nonrelativistic (Ma 1996): 



where x = k/T u , 1) , P(k, a, Q u = 0) for the pure CDM model is given by 

eq. (|I]), and d x = 0.004321, d 2 = 2.217 x lO" 6 ,^ = 11.63, and d 4 = 3.317 for k in Mpc" 1 . 
(The scale factor in T u assumes COBE normalization at a = 1. If other normalization is 
used, replace a by a<j 8 /<Ts° he .) 



Numerical simulations of structure formation in two flat C+HDM models with neutrino 
fraction Q u = 0.1 and 0.2 and two flat LCDM models with matter density fi m = 0.3 and 
0.5 are performed in order to obtain the nonlinear power spectra of matter fluctuations. 
All four simulations are performed in a (100 Mpc) 3 comoving box. The gravitational 
forces are computed with a particle-particle particle-mesh (P 3 M) algorithm (Bertschinger 
& Gelb 1991; Ma et al. 1997) with a comoving Plummer force softening length of 50 kpc. 
An identical set of random phases is used in the initial conditions for all four runs. The 
primordial power spectrum has an index of n = 1, with density fluctuations drawn from a 
random Gaussian field. A total of 128 3 simulation particles are used to represent the cold 
dark matter. For the C+HDM models, 128 3 and 10 x 128 3 particles are used to represent 
the hot component in the Q u = 0.1 and 0.2 models, respectively. Although the large particle 
number is needed to finely sample the velocity phase space (Ma & Bertschinger 1994), tests 
performed for the fl u = 0.2 model show that the power spectrum itself is little affected 




3. Nonlinear Power Spectrum 
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when 128 3 hot particles are used. Since structure forms too late in flat C+HDM models 
with Vt v > 0.2 (Ma et al. 1997 and references therein), only the Vt u = 0.1 and 0.2 models 
are studied here. Both models assume f2b = 0.05 and h = 0.5, while f2 c d m = 0.85 and 0.75 
for the two models, respectively. The two LCDM models chosen for the simulations have 
(£l m ,£l A ,h) = (0.3,0.7,0.75) and (0.5,0.5,0.7), and f2 b = 0. All four models are normalized 
to the 4-year COBE results (Bennett et al. 1996; Gorski et al. 1996). 

Figure 1 contrasts the linear (dotted) and nonlinear (solid) power spectra at various 
redshifts for the four simulated models. The hierarchical nature of gravitational collapse 
in these models is evident: the high-/c modes have become strongly nonlinear, whereas 
the low-k modes still follow the linear power spectrum. The dashed curves are from the 
analytical approximation described below. Figure 2 illustrates the dependence of the 
linear to nonlinear transformation on both time and cosmological parameters by plotting 
at various redshifts the ratio of the nonlinear and linear density variance, A n i(/c)/Ai(/c ), 
against the linear Ai(A; ) (where A = 47r/c 3 P). Note that A n i and Aj are evaluated at 
different wavenumbers, where ko — k (1 + A n i) -1 / 3 corresponds to the pre-collapsed scale of 
k, as introduced in Hamilton et al. (1991). The dependence of A nl (/c)/A 1 (A; ) on Q m , Q v , 
and time, however, demonstrates that the universality seen by Hamilton et al. is limited to 
the power-law models with n > — 1 studied there. The linear to nonlinear mapping is clearly 
more complicated in more realistic cosmological models (such as LCDM and C+HDM) 
whose spectral index changes continuously with length scale. 

Despite the apparent complication from the time dependence in Figure 2, I find 
the scaling rule given below very helpful in simplifying the construction of an analytical 
approximation. As Figure 3 shows, the time dependence is largely removed when A n i/Ai is 
plotted against the scaled quantity Ai, 

Ai=^, /3 = 0.7 + 10^ (3) 

^8 
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instead of A\. For example, (3 = 0.7 for the CDM and LCDM models, while (3 = 0.8 and 
1.1 for the Q u = 0.1 and 0.2 C+HDM models. The deviations in the two LCDM models at 
late times (a > 0.7) result from the retardation in the relative growth factor g in Q m < 1 
models. As shown in eq. (f|) below, they can be easily accounted for by including a factor 
of g 3 in the highly nonlinear regime in the fitting formula. 

The scaling behavior in eq. (j3|) can be understood in terms of the time dependence of 
the local slope of the linear power spectrum near the scale where nonlinear effects become 
important. An example is the effective spectral index n e R + 3 = din A\/dln k\k c considered 
by Jain et al. (1995), where k c is the scale at which the rms linear mass fluctuation a equals 
unity. Figure 4 illustrates the dependence of n e s on time and cosmological parameters. 
The value of n e g generally increases with time because as a grows, the wavenumber k c at 
which a = 1 decreases, and the spectral index at k c (i.e. n e s) becomes larger since the 
slope of the power spectrum for the models studied here always increases with decreasing 
k. In Figure 4, n e fr exhibits the fastest growth in the Vt v = 0.2 C+HDM model at erg > 0.3 
because the neutrino free streaming effect is more prominent in higher Q u models, which 
acts to suppress structure growth below the free streaming scale, causing Ai to bend more 
at 0.1 < k < lh Mpc^ 1 (see the dotted curves in Fig. 1). At a 8 <> 0.3 in the Vt u = 0.2 
model, on the other hand, n e g stays nearly constant because it is probing the nearly flat, 
k > 5 h Mpc -1 part of Ai. Despite this interesting behavior, Figure 4 shows that for all 
models, the dependence of n c s + 3 on as is well approximated by a power law at cr 8 0.3: 
din (n c fj + 3)/dln a 8 oc (3, where (3 = 0.7 + 10 VL 2 U as given in eq. (3). This therefore explains 
why replacing the factor n e s + 3 in the earlier work with a§ works well here. 

The simple scaling behavior introduced by eq. (|3]) allows one to approximate the 
evolution of the nonlinear power spectrum directly in terms of a§ instead of n e ff. Combining 
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these factors, I find that a close approximation for the nonlinear power spectrum is given by 



where f3 is given by eq. (3), k = k (1+A n i) 1/3 , and g = g(tt m , Q A ) and g = g(Q m (a), Uv(a)) 
are, respectively, the relative growth factor^ at present day and at a discussed in § 2. The 
time dependence is in factors o$ and g. For CDM and LCDM, a good fit is given by 
ci = 1.08 x KT 4 and c 2 = 2.10 x 10~ 5 . For C+HDM, a good fit is given by c x = 3.16 x 10~ 3 
and c 2 = 3.49 x 10~ 4 for Vt v = 0.1, and c x = 6.96 x 10~ 3 and c 2 = 4.39 x 10" 4 for Sl u = 0.2. 
The dependence of c\ and c 2 on Q u can in principle be cast in a functional form (see Ma 
1998), but since the allowed range of Q u is narrow, separate coefficients are given here in 
order to obtain the highest possible fitting accuracy. The accuracy of eq. (f|) is illustrated 
in Figure 1 (dashed curves), where the rms error for each model ranges from 3% to 10% for 
k ^ 10 h Mpc -1 at all times except z ^ 4, when the errors are about 15%. 

The functional form of G(x) in eq. (4) is chosen to give the appropriate asymptotic 
behavior A n i — > Ai in the linear regime (x<l) and A n i oc Aj^ 2 in the stable clustering 
regime (x ^> 1). In the mildly nonlinear regime, 0.1 < Ai < 1, the pre-factor [l+ln(l+0.5 x)] 
is introduced to approximate the non-negligible positive slope of A n i/Ai. This factor is 
needed because A n i and Ai are evaluated at different wavenumbers k and k , where the 
pre-collapsed ko is always smaller than k. Due to the steep positive slope of Ai in this 
region, Aj at ko is noticeably smaller than at k, and A n i(fc)/Ai(fco) is thus significantly 
above unity. Without the logarithmic pre-factor in eq. (|J) to account for this elevation, the 
approximation to A n i can be underestimated by up to 30% at 0.1 < A] < 1. 

2 Even for C+HDM models, g in eq. (4) is taken to be the familiar g(Q m (a),Q\(a)) for 
Vt v = models. The true relative growth factor for C+HDM models is in fact given by 
eq. (2), but due to its complicated scale dependence, attempts thus far to incorporate this 
factor directly in eq. (4) have not led to approximations with high accuracies (see Ma 1998). 




1 + 0.02 x 4 + ci x 8 /g 3 
1 + c 2 x 7 - 5 



(4) 



4. Summary and Discussion 



This paper presents a single formula, eq. (f|), that accurately approximates the fully 
nonlinear power spectrum of matter fluctuations for redshift z < 5 for flat CDM, LCDM, 
and C+HDM models with varying matter density Q m and neutrino fraction Q u . Eqs. (1), (2) 
and (|4]) together offer a complete description of the shape and time evolution of the matter 
power spectrum in both linear and nonlinear regime for a wide range of cosmologically 
interesting models. Figure 1 summarizes the analytical and simulation results for four 
representative models. Depending on the models and epochs, the rms errors are between 
3% and 10% for k <, 10 h Mpc" 1 at z 4 in eq. In comparison, the Peacock & Dodds 
formula has a rms error between 6% and 17% for the two LCDM simulations studied here, 
and the error reaches 50% at k ~ 10 h Mpc~ x for the two C+HDM models. 

In contrast to the scale-free models studied in Hamilton et al. (1991), the relation 
between the linear and nonlinear power spectrum is far from universal. The different panels 
of Figure 2 illustrate the dependence on cosmological parameters Q m and Q u . Moreover, 
within a given model, the fact that the curves for different times in Figure 2 do not overlap 
is important because it implies that the ratio A n i(/c)/Ai(/c ) depends not only on A\(k ) 
but also on time, or equivalently, the overall amplitude of A\(ko)- This amplitude (or time) 
dependence is present, albeit in a somewhat subtle form, in earlier work. The effective 
spectral index used by Jain et al. (1995) is clearly time-varying as shown in Figure 4 of 
this paper. The local spectral index n = d In P/d In k\k used by Peacock and Dodds (1996) 
depends on time implicitly through the factor A n i in the relation ko = k(l + A n i) -1 / 3 . 
In comparison, in this paper I have adopted the commonly used parameter og instead of 
a spectral index to characterize this time dependence, and have shown that the scaling 
behavior in eq. (|3|) and Figure 3 absorbs this dependence in C+HDM as well as LCDM 
models and results in the simple formula, eq. (4). 
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Fig. 1. — The linear and fully evolved power spectrum at various redshifts for two flat LCDM 
models (O m = 0.3 and 0.5) and two flat C+HDM models (Q v = 0.1 and 0.2). The solid curves 
are computed directly from the simulations; the dashed curves show the close approximation 
given by eq. (f|) of this paper; the dotted curves represent the linear power spectrum given 
by eqs. (1) and (2). All are normalized to the 4- year COBE data, and the present values of 
o- 8 are 1.29 (Q m = 0.3), 1.53 (Q m = 0.5), 0.9 {Q u = 0.1), and 0.81 {Q u = 0.2). In each panel, 
the curves from bottom up are for scale factors a = 0.2, 0.33, 0.6, and 1. 
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Fig. 2. — Ratio of the nonlinear to linear density variance, A nl (A;)/A 1 (A; ), as a function of 
A\(ko) at various redshifts for two LCDM and two C+HDM models. The wavenumber 
corresponds to the pre-collapsed linear value of k, where ko — k (1 + Ani)" 1 / 3 . In each panel, 
the curves from left to right correspond to a = 0.2, 0.33, 0.4, 0.6, and I. 
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Fig. 3. — Same as Figure 2, but the horizontal axis represents the scaled Ai/crf, where 
(3 = 0.7 + 10 fll- It shows that the time dependence in the ratio A n i/Ai in Figure 2 is now 
largely removed except at late times in the two LCDM models (see text). 
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Fig. 4. — The effective spectral index, n e g, as a function of as for four models. The dotted 
lines illustrate that the power law d\n(n e g + 3)/d\na 8 oc f3 is a good approximation for 
ff 8 ^ 0.3. 



